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I. INTRODUCTION 

Quantum Monte Carlo simulations (QMC) has become 
the method of choice for studying large equilibrium quan- 
tum many-body systems without approximations. While 
it is possible to obtain thermodynamic and static proper- 
ties to a high degree of accuracy with QMC, it is almost 
a paradox that estimates for the excitation spectrum 
and the equilibrium dynamics are typically obtained with 
much less accuracy. The technical reason for this is that 
QMC is invariably formulated in imaginary instead of 
real time. This is not just a matter of choice, in fact the 
imaginary time formulation is necessary to avoid crucial 
sign problems which would ruin the statistical accuracy 
of the method. The difficulty in obtaining the dynamics 
lies in transforming imaginary time correlation functions 
back to real time. This "Wick rotation" is easily carried 
out when an analytic expression of the imaginary time 
correlation function is known. However, when only nu- 
merical data and their associated error bars are available, 
as in QMC, it is well known that the direct transforma- 
tion is ill-defined and very sensitive to the errors. 

The common way to deal with this problem is to 
treat the transformation to real frequencies as a problem 
in data analysis where the imaginary time QMC result 
plays the role of the data and the real frequency spec- 
tral function is the sought-after model underlying the 
data. The data analysis problem is approached using 
Bayesian statistics which aims at identifying probabili- 
ties for different spectral functions that can account for 
the observed imaginary time data. In finding the best 
spectral function it is important that the spectral func- 
tion not only fits the data well, but also that it is con- 
sistent with prior knowledge about which types of spec- 
tral functions are permissible. The Bayesian statistical 
framework is well suited for this as both prior knowledge 
and data- fitting are taken into account. 

Although not often coined in the Bayesian language, 
the procedure of fitting certain specific functional forms 
to the imaginary time data, is an example of Bayesian 
analysis where the prior probability distribution assigns 
equal probabilities to spectral functions of the specific 
functional form and the fitting procedure selects the best 



functional parameters. However, fitting to a certain class 
of functions assumes a rather high degree of prior knowl- 
edge. While such knowledge should be used whenever 
available it is not so common that one actually knows 
the exact functional form of the spectral function a pri- 
ori. 

It is more often the case that one does not know the 
actual shape of the spectral function, but only knows 
certain sum rules and physical requirements such as real- 
valuedness and positivity. One should then prefer a prior 
probability distribution that takes only into account the 
prior knowledge and do not make extra assumptions. 
Such a maximally non-committal prior probability dis- 
tribution is gotten by maximizing the entropy of the dis- 
tribution under constraints coming from the specific a 
priori knowledgsi*^. In carrying out such a maximization 
it is important to consider the correct space to perform 
it in. A probability distribution of spectral functions is 
clearly multidimensional. Yet it is customary to treat the 
spectral function itself as a one dimensional probability 
distribution and choose a prior probability distribution 
that gives a high probability to spectral functions having 
a large entropj^. Thus instead of maximizing the en- 
tropy of the multidimensional probability distribution of 
spectral functions, the entropy of the spectral function 
itself is maximized. The latter is not the maximally non- 
committal probability distribution taking into account 
only positivity and sum-rules. In fact, to arrive at this 
socalled entropic prior involves additional assumptions^, 
which applicability to the problem at hand is question- 
able, and often one finds that methods using the entropic 
prior gives too broad spectral features. In this article we 
favor the use of another less constraining prior which re- 
flects explicitly what a priori information is included. 

In this article we use the Average Spectrum Method 
(ASM), first proposed in Ref.H, where the posterior prob- 
ability distribution is composed of a likelihood function 
and a weakly constraining prior. In the ASM the final 
spectrum is obtained as the average spectrum over this a 
posteriori probability distribution, thus the name ASM. 
We show examples of its use in getting not only the domi- 
nant features of the excitation spectra of quantum many- 
body models but also to a certain extent subdominant 
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features. 

This article is structured as follows: In section [III the 
Bayesian method is reviewed and the prior probability 
distribution is presented. The ASM is explained in Sec- 
tion [iTTI and the particular Monte Carlo implementation 
of it used in this article is described in Section HVl In sec- 
tion |V] the ASM is applied to several different quantum 
spin systems. The article ends with a summary. 

II. BAYESIAN METHOD 

The equilibrium dynamics of a physical system is char- 
acterized by the spectral function A{uj) which is real and 
non-negative. However, in QMC what is typically ob- 
tained is an imaginary time correlation function G{t) 
which is related to the spectral function as 

G(t) = J dLjK{T,Lu)A{uj) (1) 

where the kernel of the transform K{t, uj) takes on dif- 
ferent forms depending on whether the operators in the 
measured correlation function are fermionic or bosonic. 
In order to make the discussion definite and practical 
we will model the spectral function as a collection of N 
delta-functions on a frequency grid coi 

N 

AH =^A„.,5(l^-wO, (2) 

where all A^. are positive or zero. We will take a regu- 
larly spaced frequency grid such that uji is independent 
of I up to a frequency cutoff Wmax which is chosen to be 
several times the bandwidth of the system in question. 
This choice of frequency grid is not necessarily an optimal 
choice as it might be more effective to choose a finer grid 
where the spectral function is varying most. However, in 
the absence of such a priori information the choice of a 
uniform grid up to a large cutoff value is reasonable. 

Furthermore we will assume that G(r) is obtained in 
QMC simulations and recorded at discrete imaginary 
times r. With this Eq. ([1]) takes the form 

i 

The goal is to invert this relation. This is an ill-posed 
problem because of the near-zero eigenvalues of the kernel 
and therefore very sensitive to statistical errors of Gr- 
in the Bayesian approach one instead attempts to find 
the probability of a particular spectral function A given 
the QMC imaginary data G and prior knowledge. This, 
posterior probability P{A\G), can be expressed using 
Bayes theorem as 

PiA\G) cx P{G\A)P{A) (4) 

where P{G\A) is the likelihood that the QMC data turns 
out to be G given a particular spectral function A, and 



P{A) is the prior probability distribution of the spectral 
function. The prior probability distribution encodes the 
knowledge we have about the spectral function A before 
any QMC data is obtained. 

Eq. d] raises the question of how to concretely express 
the prior probability distribution P{A). We will use the 
following expression 

P(A) (X SiJ^Ko^^A^^ - Go)H,e(A^J (5) 

i 

which assigns equal probabilities to all spectral functions 
that satisfy the non-negativity requirement (Ai^. > 0) 
and the zero moment sum rule KoujiA^. = Go- In 
Eq. [n]6(a;) = 1 for a; > and zero otherwise. The prod- 
uct of 8-functions incorporates the knowledge that all 
spectral components must be non-negative, and the S- 
function constrains the spectra to obey the zero-moment 
sum rule. Higher order sum rules can be implemented 
by multiplying by more (5-functions. This prior probabil- 
ity distribution is the probability distribution having the 
highest entropy consistent with the requirement of the 
non-negativity constraint and the zeroth moment sum 
rule. It is therefore not a very selective probability dis- 
tribution as it gives the same probability to any spectral 
function that satisfy the sum rule and is non-negative. 

III. THE AVERAGE SPECTRUM METHOD 

Given the weak discriminating nature of the prior, 
Eq. ([5]), it is not a good idea to pick as the final answer 
the spectral function that maximizes the posterior proba- 
bility distribution. It is rather obvious that the spectrum 
obtained in that way will over-fit the data in the sense 
that it also will fit the noise. Instead we will pick as the 
final answer the average spectral function, obtained by 
averaging over the posterior distributioitii. Thus we will 
compute 

A==y" dAAP{A\G)/ J dAP{A\G)- (6) 

The averaging procedure itself will protect against over- 
fitting the data. The averaging procedure tends to 
smooth out the spectral function, and, in fact, it has 
been shown that when the average is carried out within 
the mean field approximation the result is identical to the 
classic MaxEnt result^. However, in general the methods 
yield different results. 

It is appropriate here to compare and contrast the 
ASM to the more commonly used MaxEnt methods^. 
The methods differ in that in MaxEnt methods an en- 
tropic prior is assumed for the spectral function and not 
the prior specified in Eq. [S] In MaxEnt methods the en- 
tropic prior is multiplied by a factor a which determines 
how much influence it has compared to the likelihood- 
function. Different MaxEnt methods differ in how the 
final answer for the spectral function is arrived at. In 
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the classic MaxEnt method the probabihty distribution 
for the parameter a, TT{a), is determined by Bayesian 
inference and the final answer is picked as the spectral 
function corresponding to the value of a that maximizes 
this probability distribution. Bryan's MaxEnt method^, 
on the other hand, is more similar to the ASM method 
as there the final spectrum is obtained by averaging the 
different spectral functions obtained at different values of 
a over 7r(a). This can either be done by computing 7r(a) 
directly for a range of a's and averaging their spectra, or 
by using a Monte Carlo procedure as shown in Ref. |8. 

Taking the average as the final answer is appropriate 
when the posterior probability has a single prominent 
peak. However, when there are more peaks the meaning 
of the average becomes more questionable. In order to de- 
tect such multiple peak situations one can focus on a few 
spectral features and make histograms of these accord- 
ing to the posterior probability distribution, and check 
for multiple peaks in these histograms. 

The averaging procedure can be efficiently carried out 
using Monte Carlo methods. In the context of getting dy- 
namics from QMC this approach is known as the Average 
Spectrum Method^, or Stochastic continuation^, but it is 
also used for data analysis in many other fields, see for 
instance Refs. [l3 and[ll|, where it is generally known as 
Markov Chain Monte Carlo methods. 

To compute the posterior probability P{A\G) we also 
need the likelihood function P{G\A). Assuming that the 
imaginary time data is distributed as Gaussians with co- 
variance matrix S, the likelihood function P{G\A) is 

P{G\A) ^ e-^T.^,(G'-G.)-i:-(G'-G.) 

where we have denoted by G' a vector of imaginary time 
values GJ. that is the average result of the i'th bin of 
QMC data containing M measurements. The assump- 
tion of having Gaussian data should be good for large 
amount of data, however this assumption should always 
be checked for instance by monitoring skewness and kur- 
tosis. Similarly we denote by Ga a vector with compo- 
nents 

j 

coming from a particular spectral function A^. In total 
there are n bins of QMC data, and for large n, S can be 
approximated by the measured covariance matrix having 
components 

"^^T. (Gr, ^.J (G;, - G,) (9) 
i 

where we have denoted by an over-bar the total mean of 
the QMC data 




It is useful to express the posterior probability in terms 
of this total mean. Using the cyclic property of the trace 
the exponent can be written as 

TrS-i^ (G' - G + G - Ga) (G' - G + G - Ga)^ 

i 

i 

+nTr (G - Ga)^ (G - Ga) . 

The first term is independent of the model A and con- 
tributes only to the normalization, thus 

PiG\A) oc e-5"Tr(G-G.)^i:-(G-G.)^ (^2) 

Note the explicit factor of n which makes the distribu- 
tion more peaked as it increases. Thus for more accurate 
QMC data (larger n) a spectral function that fits the data 
well becomes increasingly more likely than one that does 
not fit so well. This factor of n refiects the well known 
fact that the variance of the mean value is down by a 
factor 1/n. The value of n is of course rather meaning- 
less without also specifying the number of measurements 
iVnicas in each QMC bin, which determines the magni- 
tude of the components of T,. However, for a fixed large 
enough value of iVmeas, S is largely independent of n, thus 
the explicit factor of n refiects accurately how the likeli- 
hood function sharpens up when more measurements of 
QMC data is made. 

IV. MONTE CARLO IMPLEMENTATION 

The task of sampling the posterior distribution can 
be done efficiently using a Monte Carlo simulation that 
samples the distribution P{A)e~'^^^'^\ P{^) is the prior 
probability, and the energy E(A) comes from the likeli- 
hood function and is 

E{A) ^ ^fiTr {G - Ga)^ {G - Ga) , (13) 
and K = 1. 

In devising a Monte Carlo procedure one can choose 
the probability of accepting a new spectral function A' 
as 

p(A ^ A') = P(A')min(l, e-"(^(^')-^(^»). (14) 

To implement the prior probability P{A) according to 
Eq. ([5]) one starts with a spectral function that is posi- 
tive everywhere and satisfies the sum rule. In subsequent 
Monte Carlo moves one simply does not accept spectral 
functions which violate the positivity and the sum rule. 
Thus P{A) is unity for allowed spectral functions and 
zero otherwise. Typically a simulation is started with 
all spectral weight concentrated at one frequency. In 
a Monte Carlo move spectral weight is shared between 
neighboring frequencies in the following manner. First a 
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pair of neighboring frequencies uji and w^+i are chosen at 
random, and the contribution to the zero-moment sum 
rule from the spectral weights at these frequencies are 
computed: co = Ai^.Ko^^. + A^.^ji^owi+i • Then a ran- 
dom number r is selected in the interval [— co,co], and 
new spectral weights 



^U+i = A^^^,~rKoujJiKou;,+Kou,.^,) (15) 



are proposed. Note that the zero moment sum rule 
is unchanged as A^.Kq^^ + A^^^-^Kquj.+i = A'^.Kquj, + 
A'uji+i^Ouii+i- This proposed move is accepted with the 
probability specified in Eq. In particular, if either 

of the A's are negative the proposed move is rejected. 
Note that for detailed balance to hold in this scheme co 
must not change in a Monte Carlo move. For closely 
spaced frequencies this Monte Carlo move has a good ac- 
ceptance rate. To further ensure that the simulation does 
not get stuck in a local energy minimum we combine this 
move with a parallel tempering scheme in which several 
simulations of the system is simultaneously carried out at 
different temperatures l/n and a swapping move between 
different temperature configurations is included. In or- 
der to optimize the list of temperatures we have used the 
scheme in Ref. [13 where the maximum movement of con- 
figurations from the highest to the lowest temperatures 
is achieved. 

In Ref. I9I it was suggested that the entropy of the av- 
eraged spectrum be plotted vs. k and the final spectrum 
would be selected as the average at a value of k just be- 
fore the entropy makes a final drop at high values of k. 
We do not adopt such a procedure here as we find it un- 
desirable to have a procedure for selecting the spectral 
function that depends on properties of the spectral func- 
tion itself. Even though a high value of k gives solutions 
close to the most probable one there are no guarantees 
that the correct spectrum will not have a low entropy 
as is the case if the spectrum is well approximated by 
a single or a few narrow peaks. A similar criterion was 
proposed in Ref. 6 where the value of k corresponding to 
a jump in the specific heat was chosen. 

Instead we take the point of view that the final answer 
is the average spectrum at k = 1, which corresponds to 
the posterior distribution^. This means that the result- 
ing spectrum will depend on the accuracy of the input 
data, n. This is advantageous as it provides a mecha- 
nism against over-interpreting low quality data. How- 
ever, it also means that one needs to monitor how larger 
values of n will influence the flnal result. Thus a conver- 
gence analysis with n is required. This makes the method 
rather dependent on efficient QMC algorithms as gener- 
ally large values of n are needed. 



V. APPLICATIONS 

For neutron scattering the experimentally relevant 
measured quantity is the dynamic structure factor 

/oo 
dte-*(5;(i)5i,(0)) (16) 
-oc 

where the superscripts i,j indicate spin polarization di- 
rections being either x,y or 2, and S^{t) is the i'th polar- 
ization component of the spin operator in the Heisenberg 
representation at momentum q. For convenience we will 
choose units such that the lattice spacing is one. In QMC 
the accessible counterpart to the dynamic structure fac- 
tor is the imaginary time correlation function 

Sl'ir) = (5;(r)5i^(0)). (17) 

Using the Lehmann representation one finds that S*'^ and 
S'*-' are related by 

^^'^^^ ('~'^' + Sl^icj), (18) 

where /3 is the inverse temperature. Thus the kernel Kruj 
in Eq. © is 

A. Antiferromagnetic dimer in a magnetic field 

In order to check the suitability of the ASM for finding 
the spectral function we do a test on a simple system with 
a non-trivial spectrum having two peaks. We choose the 
trivial Hamiltonian of two spins in a magnetic field B 

H ^ JSi- S2- B{S^ + SI). (20) 

The dynamic structure factor of the transverse field com- 
ponents S^^{u!) displays delta-function peaks at = 
J±B each of weight tt/ [4(1 + 6' '^''{1 + 2 cosh f3B))] which 
becomes 7r/4 at low temperatures. 

We simulated this two-spin Hamiltonian at an inverse 
temperature PJ = 10 using the stochastic series expan- 
sion QMC method^ with directed loop updates^^. In the 
simulations we extracted the imaginary time correlation 
function in the x-direction at momentum vector tt. The 
imaginary time data were obtained on an equally spaced 
grid with 101 points from to (3/2, and the relative er- 
ror of the imaginary time data ranged from ^ 10~^ at 
small r to 10"^ at t — (3/2. The imaginary time data 
was then used as input to the ASM program where we 
used a regular grid with 200 frequencies having spacing 
Aw = O.OIJ. 

The results for the magnetic field value B/J = 0.1 is 
shown in Fig. [T] This result is compared to the spectrum 
obtained from the same QMC data using Bryan's Max- 
Ent method. All methods using the entropic prior gives 
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FIG. 1: Real-frequency dynamic structure factor SZ^{ijj) 
obtained from ASM (solid line) and MaxEnt (dashed line) 
for the two-spin Hamiltonian. The magnetic field value 
B/J = 0.1. The inset shows the integrated spectrum for 
the ASM curve. 



a possibility of including a default model so that the en- 
tropy is maximized when the spectral function matches 
the default model. We have used a flat model here as that 
corresponds most closely to our ASM choice of putting 
in minimal prior information. The curves in Fig. [2] were 
obtained using codes based on Ref. H. 

From Fig. [1] we see that both methods are able to re- 
solve the peaks even though the separation 2B/ J ~ 0.2. 
The peak locations corresponds well to the true value for 
both methods, but the ASM peaks are a bit narrower 
than the MaxEnt peaks. 

While the ASM gives rather sharp peaks, the two peaks 
are not equal as dictated by the exact solution. There is a 
tendency that the high energy peak is lower and broader 
than the low energy peak. This is also seen for the Max- 
Ent peaks. The spectral weight is however equally dis- 
tributed on the two peaks in both the low and the high 
field cases, see inset of Fig. [TJ We expect that the peaks 
become more and more equal as the quality of the QMC 
data is increased (larger n). This has the effect that the 
likelihood function becomes more peaked and more de- 
tails of the spectrum will be better resolved. An example 
of this is shown in Fig. [2] where it is clear the the dou- 
ble peak structure is only revealed for data of sufficient 
quality. 

We have also simulated the dimer system with a bigger 
value of the magnetic field, B=0.4J. For this value of B 
the peaks aX to — 0.6 J and lj — lAJ are very narrow in 
both the ASM and the Bryan MaxEnt method. 
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FIG. 2: (Color online) The effect of improving the data qual- 
ity by increasing the number of Monte Carlo bins n. Each 
panel shows the dynamic structure factor for B = O.IJ for 
three independent data set (different line styles). The num- 
ber of data bins were n = 4 (left), n = 10 (middle), and 
n = 20 (right). For comparison the results shown in Fig. [T] 
was carried out using n — 200. 



B. Spin-1 chain 

We now move on to a nontrivial example, the spin- 
1 antiferromagnetic chain, the so called Haldane chain. 
The Haldane chain is famous for being gapped in con- 
trast to the half-integer spin chains^^. The minimum 
gap is at Q = TT in units of the inverse lattice spacing. 
Fig. [3] shows plots of SqL^ (uj) for different temperatures 
obtained using the ASM. Note how the peak position 
and width increase with temperature. To compare with 
MaxEnt we have shown the MaxEnt result using Bryan's 
method for a single temperature T/J — 0.25 as a dashed 
curve. Note that the MaxEnt curve captures the peak po- 
sition well, but gives a very broad peak. The inset shows 
a comparison of the temperature dependence of the gap 
vs. a non-linear sigma model prediction which was ob- 
tained by solving the finite temperature gap equation in 
Ref. [161 numerically. In the inset we also show a compar- 
ison of the width of the peaks, quantified by their full 
width at half maximum (FWHM), with predicted val- 
ues from a combined nonlinear cr-model and scattering 
matrix calculationii. The agreement is quite remarkable 
and involves no adjustable parameters. 

One can ask whether the temperature broadening of 
the peak seen in Fig. [3] obtained using the ASM is just 
due to the "motion" of a single sharp peak. Fig. 0] shows 
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FIG. 3: (Color online) Dynamic structure factor for 
a ID spin chain with 64 sites obtained from the ASM (solid 
lines) at different temperatures indicated by the legends. The 
dashed curve is the MaxEnt result for T/ J — 0.25. The curves 
for T/J = 0.0625 and T/J = 0.125 have been scaled down 
by a factor 1/2 to fit inside the figure boundaries. The inset 
shows the peak positions A (circles) and peak widths FWHM 
(triangles) as functions of temperature. The solid lines are 
the (T-model predictions for these quantities. 
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FIG. 4: Snapshots of spectra. These spectra (and others) 
are averaged over in order to yield the result shown in Fig. [3] 
The spectra here are all for T/J — 0.25. 



that this is not the case. 



C. Bond alternating antiferromagnetic chain 

Another nontrivial spin model is the bond alternat- 
ing spin-1/2 Heisenberg chain (BAHC) which has been 
studied extensively, and is relevant for materials such 
as Cu(N03)2 • 2.5D2Oi^>i^'20, the spin-Peirels material 
CuGeOa^ and others (See Refi^^). The Hamiltonian for 




Q [7i/a] 



FIG. 5: Gray scale plot of Sq (lu) for the BAHC with A = 0.8. 
The simulations were carried out at I3J = 16 on a lattice with 
128 sites and periodic boundary conditions. The solid blue 
curve indicates the one-magnon excitations as calculated us- 
ing a series expansion about the dimer limit^, and the dotted 
lines show the kinematic boundaries of two-particle excita- 
tions. 



the BAHC is 

H = JJ2 (^2^-1 • 52. + XS2^ ■ S2^+l) (21) 

i 

where A > 0. Although the BAHC is a one-dimensional 
model, it is not solvable by the Bethe Ansatz. Thus 
other techniques are needed to obtain the dynamics. In 
this regard investigations using bosonization^^ the RPA 
approximationSi, series expansion a^^'^^i^^'^^i^^'^^i'^° and 
exact diagonalization studies'^'' have produced very im- 
pressive results for the dynamics of the BAHC containing 
predictions of the dispersion of one magnon excitations 
as well as bound states and details about multi particle 
excitations. 

We carried out QMC simulations of the BAHC for a 
chain with 128 sites and periodic boundary conditions 
at inverse temperature (3J = 16 and A = 0.8. The 
ASM was used to obtain the spectra at all momentum 
points. Figure [5] shows a gray scale plot of Sq{uj) for 
different values of Q and lu . The one magnon excita- 
tions are easily identified as the sharp dark feature and 
agrees very well with that obtained from series expan- 
sion to order A^^^ shown as the blue solid curve. For 
Q ^ 0.57r many-particle excitations are visible. This 
agrees qualitatively with the results in Ref. 30 which 
shows that the many-particle continuum has apprecia- 
bly more spectral weight for Q > G.Stt than for smaller 
Q. For O.Stt ^ Q ^ O.TStt there is an almost flat fea- 
ture in the continuum at a; ~ 1.9J which is well sepa- 
rated from the band of one magnon excitations and also 
from the kinematic boundaries of two magnon excitations 
shown as blue dotted lines. This is not seen from the se- 



ries expansion^° and RPA results which predict that 
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FIG. 6: Line shapes at a fixed momentum Q — Sn/A for QMC 
data sets of different lengtlis n indicated by tlie legends. 
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FIG. 7; Gray scale plot of Sq" (lj) for the BAHC with A = 0.8 
in a magnetic field B — 0.2J. The inverse temperature is 
l3J — 16 and L — 128. The solid lines are the spin-split one 
magnon result. 



the continuum should have biggest spectral weight at its 
lower boundary. However, this feature is reminiscent of 
that seen in experiments on Cu(N03)2 -2.5020^ where a 
dispersion-less feature in the continuum was reported. As 
Q is increased towards tt this feature broadens and van- 
ishes. Some structure reappears in the continuum close 
to Q TT where a peak at lu ^ J and a very weak feature 
at a; ^ 2 J is seen. A word of caution is needed in inter- 
preting weak features of Fig.O This is because Fig.[5]also 
shows occurrence of spectral weight in between the one 
magnon peak and the lower kinematic boundary of the 
two magnon excitations, where one expects a gap. This 
is probably caused by insufficient quality of the QMC 
data which gives spectral weight in unwanted places in 
a similar fashion to what is seen in Fig. [1] at a; J for 
B = O.IJ. 

The QMC data plotted in Fig. [5] were taken from a 
run with in all n — 2000 data bins. In order to see how 
the number of QMC bins affect the line shapes we show 
in Fig. [6] the dynamic structure factor at Q = 3tt/A for 
three different values of n. While there is some signifi- 
cant change in the line shape from n = 20 to n — 200, 
increasing n to 2000 has only minor effects. 

We will now add a magnetic field term —Bj^i Sf to 
Eq. dH])- For A = the BAHC is just a collection of 
independent antiferromagnetic dimers. When subjecting 
a dimer to a magnetic field in the spin z direction the 
degeneracy of the spin triplet excitations is lifted, and 
one expects a double-peak structure, as seen in Fig.[Tl in 
the transverse dynamic structure factor S^^ . For finite A 
the dimers become coupled, however one still expects the 
splitting to occur, at least for small values of the mag- 
netic field. Fig. [7] shows a gray scale plot of S''q{lo) for 
A = 0.8 and a small value of the magnetic field B = 0.2 J. 



The splitting of the one magnon peak is clearly seen and 
agrees, for small values of Q, well with the expectation 
that the effect of the magnetic field is simply to displace 
the one magnon dispersion by ±B. The solid lines indi- 
cate this. We have taken the one magnon dispersion from 
the series expansion^^ and added (subtracted) an energy 
B = 0.2J. For O.Stt < Q < 0.757r there are deviations 
from this simple picture, as the upper branch is higher 
in energy and broadens considerably. For even higher 
momentum values there is significant broadening of the 
peaks and at Q = tt they are hardly distinguishable. For 
Q ^ 0.757r one can also see the appearance of many- 
particle excitations above the one magnon peaks. 

For a large value of the magnetic field the lower branch 
goes to zero energy at a certain characteristic value of the 
momentum. Figure [S] shows a gray scale plot of the trans- 
verse structure factor Sq^{uj) for A = 0.8 and B = J. 
One can clearly see that there is a branch of excitations 
that approaches zero at Q « O.Stt and at Q = tt. This 
is consistent with the results reported in Ref. [sl. It is 
also apparent that the intensity at Q ~ O.Stt vanishes as 
the energy approaches zero, while the intensity at Q — tt 
is high. The high energy magnon branch is clearly seen 
for Q < O.Btt and gets broadened considerably and dis- 
appears for larger Q. There is also a sharp finite energy 
peak seen at small Q resulting from the merger of the 
two magnon branches. 



D. Heisenberg antiferromagnetic chain 

The spin-1/2 Heisenberg chain was the first nontrivial 
quantum many-body problem to be solved exactly^S,. Yet 
it is still only recently that exact results for the dynam- 
ical correlation functions have appeareci^. We compare 
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FIG. 8: Gray scale plot of S"''=(Q,w) for the BAHC in a 
magnetic field B = J. A = 0.8, I3J = 16 and L = 128. 

here the ASM with the exact numerical result for the dy- 
namic structure factor for the Heisenberg antiferromag- 
netic chain. 

In Fig. [9] we show the lineshape of S'^^(Q, at Q = 
O.Stt, where the gap is the largest, as well as at Q = O.Ott 
where the exact result has a very long high-energy tail. 
We see that the exact results (red dashed curves) is zero 
up to a certain energy where a vertical leading edge marks 
the onset of a continuum of excitations. The ASM results 
have no true vertical leading edge, but rather a power- 
law increase. This smooth increase is inevitable in the 
ASM method as even a prior that incorporates a strict 
requirement of having a vertical leading edge will give 
a smooth leading edge if there is uncertainity about the 
position of the edge. There is also a slight difference in 
the location of the maximum intensity. While the exact 
results peak right at the leading edge, the ASM results 
peak slightly above the exact results. This is most promi- 
nent in the Q = O.Stt case and is probably because the 
true lineshapes are very asymmetric and tend to push up 
the peak in energy. This asymmetry can also be seen in 
both ASM curves. The extent of how high up in energy 
the continuum reaches can be seen from the insets. The 
high energy tail is very well reproduced by the ASM for 
Q = O.Ott while it is overestimated for the Q = O.Stt case. 



E. Square lattice Heisenberg antiferromagnet 

The spin-1/2 square lattice Heisenberg antiferromag- 
net (2DAF) has been studied intensively because of its 
relevance to the cuprate materials that are superconduct- 
ing at high temperatures when doped. The dynamics of 
the 2DAF is rather well described by linear spin-wave 
theory'^^. However, linear spin wave theory does not ac- 
count for a magnon dispersion along the zone bound- 




FIG. 9: (Color online) Line shapes of S''(Q,uj) for the ID 
Heisenberg antiferromagnet at Q = O.Stt (upper panel) and 
Q = O.Qtf (lower panel) solid curves. The red dashed lines 
are exact results obtained from the Bethe Ansatz. The chain 
has periodic boundary conditions and has L = 500 sites. The 
QMC simulations are carried out at j3J — 40 while the Bethe 
Ansatz result is obtained at T = 0. The insets shows the 
same results but on a semi-log scale. 



ary. Such a dispersion was predicted using an expan- 
sion around the Ising limi t^^'^^ and indicates a differ- 
ence in energy between the magnon peaks at (tt, 0) and 
(tt/2, Tr/2) of about 7-9%, the energy at (tt/2, tt/2) being 
the highest. Similar result was obtained using QMC: In 
Ref. '37 the QMC data were fitted to a functional form 
consisting of a delta-function and a broad continuum, 
while in Ref. [H the MaxEnt method was used. Higher 
order Holstein-Primakoff spin wave calculations gives a 
smaller value, 2%^^, as does an expansion based on the 
Dyson-Maleev transformation^Siii. 

Experimental measurements of the material copper 
formate tetradeuterate(CFTD)^'^-'^*? indicated a differ- 
ence of 7% in agreement with the series expansion re- 
sults and the QMC, however La2Cu04 shows^ an en- 
tirely different dispersion with the peak at (tt, 0) being 
higher in energy than at (7t/2,tt/2). This dispersion 
has been explained as special features of the Hubbard 
model^. Recently experiments on K2V3O8, also sup- 
posedly a realization of the Heisenberg antiferromagnet 
on the square lattice, showed a double peak structure 
of unknown origin at (tt/2,tt/2)'^^. In order to inves- 
tigate this possible double peak structure we repeated 
the simulations of Ref. ^33 and analyzed the imaginary 
time data using the ASM which gives unbiased infor- 
mation about the line shapes. In order to distinguish 
transversal and longitudinal excitations the simulations 
were carried out as in Ref. [13 by imposing a staggered 
magnetic field -ffstag = 0.00161S that yields a staggered 
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FIG. 10: Transverse (solid curves) and longitudinal (dashed 
curves) dynamic structure factor for the 2DAF at Q = (tt, 0) 
(upper panel) and Q — (7r/2,7r/2) (lower panel). 



magnetization consistent with the experimental value 
nis = 0.307 on a 32 X 32 lattice at an inverse temper- 
ature f3J = 32. We measured both the transverse dy- 
namic structure factor S^^ and the longitudinal one S^^ . 
The results for the two momentum points Q = (tt, 0) 
and Q — (7r/2,7r/2) are shown in Fig. [TOl We observe a 
difference in magnon energies in the transverse channel 



corresponding to (£;(^/2,^/2) - £^(7r,o))/^^(7r 



/2,7r/2) ~ 6%, 

determined from the location of the maximum. However, 
the peak locations are at slightly higher energies than the 
corresponding delta- function locations found in Ref. [s^- 
As we expect a priori that the dynamic structure factor in 
the transverse channel contains a delta- function like one 
magnoii peak and a continuum we believe that the result 
in Ref. |37| is the most accurate as it accounts for more 
prior information. However, for the longitudinal chan- 
nel the expected functional form of the spectral function 
is not so clear. In particular it is not obvious that the 
particular functional form chosen in Ref. [s^ in the lon- 
gitudinal channel is flexible enough to track the real line 
shape. In fact, in contrast to the result reported there, 
at Q = (7r/2,7r/2), the lower panel of Fig. [TOl shows that 
the peak location in the transverse channel is at a sub- 
stantial lower energy {^-^ 10%) than the peak in the longi- 
tudinal channel. For an experiment that measures both 
the longitudinal and transverse structure factors simulta- 
neously this could give rise to a double peak structure at 
(7r/2, 7r/2). Such a double peak should also be apparent 
at (tTjO), although more weakly, because the longitudi- 
nal structure factor is more strongly peaked at (7r/2, 7r/2) 
than at (7r,0). In fact, as can be seen from Fig. [TUl the 
longitudinal dynamic structure factor at (tt, 0) has a very 
long high-energy tail. 



VI. SUMMARY 

Obtaining equilibrium dynamics from numerical imag- 
inary time correlation functions is an important task. We 
have in this article investigated the suitability of a specific 
Bayesian method for doing this. This method, known as 
the ASM^, proposed already in 1991 has not been widely 
used. We suspect that this is because its nature is such 
that for it to give good results one needs rather accu- 
rate QMC data. However, QMC simulations have im- 
proved considerably the last years, thus it is timely to 
reconsider its usefulness. The ASM is a Bayesian data 
analysis method, where instead of picking the final result 
as the spectrum that maximizes the posterior probabil- 
ity distribution, the final answer is picked as the averaged 
spectrum over the posterior probability distribution. The 
reason for selecting to take the average is the rather un- 
selective nature of the specific prior probability distribu- 
tion used. We argue for the use of a prior probability 
distribution that encodes just hard knowledge; spectral 
positivity and sum rules, and the specific form of the 
prior is then the one maximizing the information theory 
entropy under these constraints. One should note that 
this prior is not the entropic prior used in various Max- 
Ent methods. The entropic prior gives high probabilities 
to spectral functions that itself has high entropy, thus 
favoring smooth spectral functions. 

There are other methods that resembles the ASM. The 
Stochastic continuation method^ is essentially the same 
method, except for the use of a drop in entropy as the 
criterion for determining the temperature at which the 
sampling is carried out. In the ASM the posterior prob- 
ability distribution is sampled directly. Thus in essence 
the quality of the input data determines the effective sam- 
pling temperature which is implicit in the approach. We 
find this desirable as it protects from over interpreting 
bad data and makes the procedure independent of the 
particular form of the spectral function itself. However, 
this also implies the need of a convergence analysis of the 
obtained spectral function with increasingly better QMC 
data. Some MaxEnt methods, such as the Bryan Max- 
Ent method, also outputs as the final answer an averaged 
spectrum. In the case of Bryan's method^ the average is 
taken over the probability distribution of the coefficient 
determining the relative importance of the entropic prior. 

The ASM is on at least as firm statistical footings as 
other Bayesian methods'^. It has the disadvantage of 
being computationally demanding, however it is not as 
computer-intensive as the QMC simulations themselves. 
A typical run of the ASM, for one momentum space 
point, takes about 4 hours on an Intel Pentium IV, 2.4 
GHz processor. In comparison running MaxEnt methods 
takes typically of the order of tens of seconds. 

In showing examples of the ASM we have sampled the 
posterior probability distribution and obtained spectral 
functions for several model systems. Of new results we 
have shown that using this method we can obtain the 
finite temperature position and broadening of the Hal- 
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dane gap in spin-1 antiferromagnets, and that the results 
agree very well with nonlinear cr-model predictions with- 
out any adjustable parameters. We have also applied the 
method to the spin- 1/2 Heisenberg chain with alternat- 
ing bond strengths where we found a quantitative very 
good agreement with other methods for the dispersion of 
one magnon excitations. We also observed some struc- 
ture in the continuum of many-particle excitations which 
have not been seen using other methods. At present it 
is unclear whether these many-particle features are real 
or whether they are artifacts of insufficient QMC data. 
We have also added a magnetic field to the bond alter- 
nating chain, and observed the expected spin-split spec- 
trum in the transverse dynamic structure factor. For a 
bigger value of the magnetic field we also see the weak 
incommensurate low-energy mode and the much stronger 
low-energy mode ai Q — it. We have compared the ASM 
for the ID Heisenberg antiferromagnet with exact Bethe 
Ansatz results. The comparison reveals a good similarity, 
although certain sharp features of the exact result, such 
as the lineshape's vertical leading edge, is not accurately 
reproduced by the ASM. Finally we studied the dynamic 



structure factor at the zone boundary for the two dimen- 
sional square lattice spin- 1/2 Heisenberg antiferromag- 
net, and found results consistent with existing results on 
that system except for a difference in peak locations of 
the transverse and longitudinal dynamic structure fac- 
tor at the same momentum value that can possibly give 
rise to a double peak structure in measurements using 
unpolarized neutrons. 
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